module GOLD_PressureForce_AFV
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of GOLD.                                        *
!*                                                                     *
!* GOLD is free software; you can redistribute it and/or modify it and *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* GOLD is distributed in the hope that it will be useful, but WITHOUT *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg and Alistair Adcroft, 2006                      *
!*                                                                     *
!*    This file contains the subroutine that calculates the hori-      *
!*  zontal accelerations due to pressure gradients using a 2nd-order   *
!*  analitically vertically integrated finite volume form, as          *
!*  described by Adcroft, A., R. Hallberg, and M. Harrison, 2008: A    *
!*  finite volume discretization of the pressure gradient force using  *
!*  analytic integration. Ocean Modelling, 22, 106-113.  This form     *
!*  eliminates the thermobaric instabilities that had been a problem   *
!*  with previous forms of the pressure gradient force calculation,    *
!*  as described by Hallberg (2005), Ocean Modelling, 8, 279-300.      *
!*                                                                     *
!*    PressureForce takes 8 arguments, which are described below. If   *
!*  a non-split time stepping scheme is used, the last three arguments *
!*  are ignored.                                                       *
!*                                                                     *
!*  Macros written all in capital letters are defined in GOLD_memory.h.*
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, f                                     *
!*    j+1  > o > o >   At ^:  v, PFv                                   *
!*    j    x ^ x ^ x   At >:  u, PFu                                   *
!*    j    > o > o >   At o:  h, D, M, e, p, pbce, gtot, T, S, rho_st  *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1                                                  *
!*           i  i+1                                                    *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use GOLD_diag_mediator, only : post_data, register_diag_field
use GOLD_diag_mediator, only : safe_alloc_ptr, diag_ptrs, time_type
use GOLD_error_handler, only : GOLD_error, FATAL, WARNING, is_root_pe
use GOLD_file_parser, only : get_param, log_param, log_version, param_file_type
use GOLD_grid, only : ocean_grid_type
use GOLD_PressureForce_Mont, only : set_pbce_Bouss, set_pbce_nonBouss
use GOLD_tidal_forcing, only : calc_tidal_forcing, tidal_forcing_CS
use GOLD_variables, only : thermo_var_ptrs
use GOLD_EOS, only : calculate_density, calculate_density_derivs
use GOLD_EOS, only : int_density_dz, int_specific_vol_dp
implicit none ; private

#include <GOLD_memory.h>

public PressureForce_AFV, PressureForce_AFV_init, PressureForce_AFV_end
public PressureForce_AFV_Bouss, PressureForce_AFV_nonBouss

type, public :: PressureForce_AFV_CS ; private
  logical :: tides          ! If true, apply tidal momentum forcing.
  real    :: Rho0           !   The density used in the Boussinesq
                            ! approximation, in kg m-3.
  real    :: GFS_scale      !   A scaling of the surface pressure gradients to
                            ! allow the use of a reduced gravity model.
  type(time_type), pointer :: Time ! A pointer to the ocean model's clock.
  type(diag_ptrs), pointer :: diag ! A pointer to a structure of shareable
                             ! ocean diagnostic fields.
  integer :: id_e_tidal = -1
  type(tidal_forcing_CS), pointer :: tides_CSp => NULL()
end type PressureForce_AFV_CS

contains

subroutine PressureForce_AFV(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
  real, dimension(NXMEM_,NYMEM_,NZ_), intent(in)   :: h
  type(thermo_var_ptrs), intent(inout)             :: tv
  real, dimension(NXMEMQ_,NYMEM_,NZ_), intent(out) :: PFu
  real, dimension(NXMEM_,NYMEMQ_,NZ_), intent(out) :: PFv
  type(ocean_grid_type),               intent(in)  :: G
  type(PressureForce_AFV_CS),              pointer     :: CS
  real, dimension(:,:),               optional, pointer     :: p_atm
  real, dimension(NXMEM_,NYMEM_,NZ_), optional, intent(out) :: pbce
  real, dimension(NXMEM_,NYMEM_),     optional, intent(out) :: eta

!    This subroutine works as a temporary interface between the model and the
! Boussinesq and non-Boussinesq pressure force routines.
! Descriptions of the variables are in each of the routines called in the
! following conditional block.

  if (G%Boussinesq) then
    call PressureForce_AFV_bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
  else
    call PressureForce_AFV_nonbouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
  endif

end subroutine PressureForce_AFV

subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(in)    :: h
  type(thermo_var_ptrs),               intent(inout) :: tv
  real, dimension(NXMEMQ_,NYMEM_,NZ_), intent(out)  :: PFu
  real, dimension(NXMEM_,NYMEMQ_,NZ_), intent(out)  :: PFv
  type(ocean_grid_type),               intent(in)   :: G
  type(PressureForce_AFV_CS),          pointer      :: CS
  real, dimension(:,:),               optional, pointer     :: p_atm
  real, dimension(NXMEM_,NYMEM_,NZ_), optional, intent(out) :: pbce
  real, dimension(NXMEM_,NYMEM_),     optional, intent(out) :: eta

!    This subroutine determines the acceleration due to hydrostatic pressure
!  forces, using the analytic finite volume form of the Pressure gradient, and
!  does not make the Boussinesq approximation. It follows the approach described
!  in Adcroft, Hallberg, and Harrison, 2008, Ocean Modelling.
!    To work, the following fields must be set outside of the usual
!  ie to ie, je to je range before this subroutine is called:
!   h[ie+1] and h[je+1] and (if tv%eqn_of_state is set) T[ie+1], S[ie+1],
!   T[je+1], and S[je+1].
! Arguments: h - Layer thickness, in H (probably kg m-2).
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (out)     PFu - Zonal acceleration due to pressure gradients
!                  (equal to -dM/dx) in m s-2.
!  (out)     PFv - Meridional acceleration due to pressure
!                  gradients (equal to -dM/dy) in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 PressureForce_init.
!  (in)      p_atm - The pressure at the ice-ocean or atmosphere-ocean
!                    interface in Pa.
!  (out)     pbce - The baroclinic pressure anomaly in each layer
!                   due to eta anomalies, in m2 s-2 H-1.
!                   pbce points to a space with nz layers or NULL.
!  (out)     eta - The bottom mass used to calculate PFu and PFv, in H,
!                  with any tidal contributions or compressibility compensation.

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p ! Interface pressure in Pa.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
    T_tmp, &    ! Temporary array of temperatures where layers that are lighter
                ! than the mixed layer have the mixed layer's properties, in C.
    S_tmp       ! Temporary array of salinities where layers that are lighter
                ! than the mixed layer have the mixed layer's properties, in psu.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: &
    dza, &      ! The change in geopotential anomaly between the top and bottom
                ! of a layer, in m2 s-2.
    intp_dza    ! The vertical integral in depth of the pressure anomaly less
                ! the pressure anomaly at the top of the layer, in Pa m2 s-2.
  real, dimension(SZI_(G),SZJ_(G))  :: &
    dp, &       ! The (positive) change in pressure across a layer, in Pa.
    SSH, &      ! The sea surface height anomaly, in m.
    e_tidal, &  ! The bottom geopotential anomaly due to tidal forces from
                ! astronomical sources and self-attraction and loading, in m.
    dM, &       ! The barotropic adjustment to the Montgomery potential to
                ! account for a reduced gravity model, in m2 s-2.
    Rho_cv_BL,& !  The coordinate potential density in the deepest variable
                ! density near-surface layer, in kg m-3.
    za          ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the
                ! interface atop a layer, in m2 s-2.
  real, dimension(SZIQ_(G),SZJ_(G)) :: &
    intx_za     ! The zonal integral of the geopotential anomaly along the
                ! interface below a layer, divided by the grid spacing, m2 s-2.
  real, dimension(SZIQ_(G),SZJ_(G),SZK_(G)) :: &
    intx_dza    ! The change in intx_za through a layer, in m2 s-2.
  real, dimension(SZI_(G),SZJQ_(G)) :: &
    inty_za     ! The meridional integral of the geopotential anomaly along the
                ! interface below a layer, divided by the grid spacing, m2 s-2.
  real, dimension(SZI_(G),SZJQ_(G),SZK_(G)) :: &
    inty_dza    ! The change in inty_za through a layer, in m2 s-2.
  real :: p_ref(SZI_(G))     !   The pressure used to calculate the coordinate
                             ! density, in Pa (usually 2e7 Pa = 2000 dbar).

  real :: dp_neglect         ! A thickness that is so small it is usually lost
                             ! in roundoff and can be neglected, in Pa.
  real :: alpha_anom         ! The in-situ specific volume, averaged over a
                             ! layer, less alpha_ref, in m3 kg-1.
  real :: aaL, aaR           ! The in-situ specific volume at the thickness point
                             ! to the" left" or "right" of a velocity point,
                             ! averaged over a layer, less Rho_ref, in kg m-3.
  logical :: use_p_atm       ! If true, use the atmospheric pressure.
  logical :: use_EOS    ! If true, density is calculated from T & S using an
                        ! equation of state.
  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.

  real :: alpha_ref ! A reference specific volume, in m3 kg-1, that is used
                    ! to reduce the impact of truncation errors.
  real :: rho_in_situ(SZI_(G)) ! The in situ density, in kg m-3.
  real :: Pa_to_H   ! A factor to convert from Pa to the thicknesss units (H).
!  real :: oneatm = 101325.0  ! 1 atm in Pa (kg/ms2)
  real :: I_gEarth
  real, parameter :: C1_6 = 1.0/6.0
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
  integer :: i, j, k
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  nkmb=G%nk_rho_varies
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq

  use_p_atm = .false.
  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
  use_EOS = associated(tv%eqn_of_state)

  if (.not.associated(CS)) call GOLD_error(FATAL, &
       "GOLD_PressureForce: Module must be initialized before it is used.")

  dp_neglect = G%H_to_Pa * G%H_subroundoff
  alpha_ref = 1.0/CS%Rho0

  if (use_p_atm) then
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      p(i,j,1) = p_atm(i,j)
    enddo ; enddo
  else
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      p(i,j,1) = 0.0 ! or oneatm
    enddo ; enddo
  endif
  do k=2,nz+1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
    p(i,j,K) = p(i,j,K-1) + G%H_to_Pa * h(i,j,k-1)
  enddo ; enddo ; enddo

  I_gEarth = 1.0 / G%g_Earth

  if (use_EOS) then
  !   With a bulk mixed layer, replace the T & S of any layers that are
  ! lighter than the the buffer layer with the properties of the buffer
  ! layer.  These layers will be massless anyway, and it avoids any
  ! formal calculations with hydrostatically unstable profiles.
    if (nkmb>0) then
      tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp
      tv_tmp%eqn_of_state => tv%eqn_of_state
      do k=1,nkmb ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
      enddo ; enddo ; enddo
      do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
      do j=Jsq,Jeq+1
        call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
                        Rho_cv_BL(:,j), Isq, Ieq-Isq+2, tv%eqn_of_state)
      enddo
      do k=nkmb+1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        if (G%Rlay(k) < Rho_cv_BL(i,j)) then
          tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
        else
          tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
        endif
      enddo ; enddo ; enddo
    else
      tv_tmp%T => tv%T ; tv_tmp%S => tv%S
      tv_tmp%eqn_of_state => tv%eqn_of_state
    endif
  endif

  do k=1,nz
    ! Calculate 4 integrals through the layer that are required in the
    ! subsequent calculation.
    if (use_EOS) then
      call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,K), &
                               p(:,:,K+1), alpha_ref, G, tv%eqn_of_state, &
                               dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
                               inty_dza(:,:,k))
    else
      alpha_anom = 1.0/G%Rlay(k) - alpha_ref
      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        dp(i,j) = G%H_to_Pa * h(i,j,k)
        dza(i,j,k) = alpha_anom * dp(i,j)
        intp_dza(i,j,k) = 0.5 * alpha_anom * dp(i,j)**2
      enddo ; enddo
      do j=js,je ; do I=Isq,Ieq
        intx_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i+1,j))
      enddo ; enddo
      do J=Jsq,Jeq ; do i=is,ie
        inty_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i,j+1))
      enddo ; enddo
    endif
  enddo

  !   The bottom geopotential anomaly is calculated first so that the increments
  ! to the geopotential anomalies can be reused.  Alternately, the surface
  ! geopotential could be calculated directly with separate calls to
  ! int_specific_vol_dp with alpha_ref=0, and the anomalies used going
  ! downward, which would relieve the need for dza, intp_dza, intx_dza, and
  ! inty_dza to be 3-D arrays.

  ! Sum vertically to determine the surface geopotential anomaly.
  do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
    za(i,j) = alpha_ref*p(i,j,nz+1) - G%g_Earth*G%D(i,j)
  enddo ; enddo
  do k=nz,1,-1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
    za(i,j) = za(i,j) + dza(i,j,k)
  enddo ; enddo ; enddo

  if (CS%tides) then
    ! Find and add the tidal geopotential anomaly.
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth
    enddo ; enddo
    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp)
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      za(i,j) = za(i,j) - G%g_Earth*e_tidal(i,j)
    enddo ; enddo
  endif

  if (CS%GFS_scale < 1.0) then
    ! Adjust the Montgomery potential to make this a reduced gravity model.
    if (use_EOS) then
      do j=Jsq,Jeq+1
        call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), &
                               rho_in_situ, Isq, Ieq-Isq+2, tv%eqn_of_state)

        do i=Isq,Ieq+1
          dM(i,j) = (CS%GFS_scale - 1.0) * &
            (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
        enddo
      enddo
    else
      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        dM(i,j) = (CS%GFS_scale - 1.0) * &
          (p(i,j,1)*(1.0/G%Rlay(1) - alpha_ref) + za(i,j))
      enddo ; enddo
    endif
!  else
!    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; dM(i,j) = 0.0 ; enddo ; enddo
  endif

  !   This order of integrating upward and then downward again is necessary with
  ! a nonlinear equation of state, so that the surface geopotentials will go
  ! linearly between the values at thickness points, but the bottom
  ! geopotentials will not now be linear at the sub-grid-scale.  Doing this
  ! ensures no motion with flat isopycnals, even with a nonlinear equation of state.
  do j=js,je ; do I=Isq,Ieq
    intx_za(I,j) = 0.5*(za(i,j) + za(i+1,j))
  enddo ; enddo
  do J=Jsq,Jeq ; do i=is,ie
    inty_za(i,J) = 0.5*(za(i,j) + za(i,j+1))
  enddo ; enddo
  do k=1,nz
    ! These expressions for the acceleration have been carefully checked in
    ! a set of idealized cases, and should be bug-free.
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      dp(i,j) = G%H_to_Pa*h(i,j,k)
      za(i,j) = za(i,j) - dza(i,j,k)
    enddo ; enddo
    do j=js,je ; do I=Isq,Ieq
      intx_za(I,j) = intx_za(I,j) - intx_dza(I,j,k)
      PFu(I,j,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - &
                     (za(i+1,j)*dp(i+1,j) + intp_dza(i+1,j,k))) + &
                    ((dp(i+1,j) - dp(i,j)) * intx_za(I,j) - &
                     (p(i+1,j,K) - p(i,j,K)) * intx_dza(I,j,k))) * &
                   (2.0*G%IDXu(I,j) / ((dp(i,j) + dp(i+1,j)) + dp_neglect))
    enddo ; enddo
    do J=Jsq,Jeq ; do i=is,ie
      inty_za(i,J) = inty_za(i,J) - inty_dza(i,J,k)
      PFv(i,J,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - &
                     (za(i,j+1)*dp(i,j+1) + intp_dza(i,j+1,k))) + &
                    ((dp(i,j+1) - dp(i,j)) * inty_za(i,J) - &
                     (p(i,j+1,K) - p(i,j,K)) * inty_dza(i,J,k))) * &
                   (2.0*G%IDYv(i,J) / ((dp(i,j) + dp(i,j+1)) + dp_neglect))
    enddo ; enddo

    if (CS%GFS_scale < 1.0) then
      ! Adjust the Montgomery potential to make this a reduced gravity model.
      do j=js,je ; do I=Isq,Ieq
        PFu(I,j,k) = PFu(I,j,k) - (dM(i+1,j) - dM(i,j)) * G%IDXu(I,j)
      enddo ; enddo
      do J=Jsq,Jeq ; do i=is,ie
        PFv(i,J,k) = PFv(i,J,k) - (dM(i,j+1) - dM(i,j)) * G%IDYv(i,J)
      enddo ; enddo
    endif
  enddo

  if (present(pbce)) then
    call set_pbce_nonBouss(p, tv_tmp, G, G%g_Earth, CS%GFS_scale, pbce)
  endif

  if (present(eta)) then
    Pa_to_H = 1.0 / G%H_to_Pa
    if (use_p_atm) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*Pa_to_H ! eta has the same units as h.
    enddo ; enddo ; else ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      eta(i,j) = p(i,j,nz+1)*Pa_to_H ! eta has the same units as h.
    enddo ; enddo ; endif
  endif

  if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag)

end subroutine PressureForce_AFV_nonBouss

subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(in)    :: h
  type(thermo_var_ptrs),               intent(inout) :: tv
  real, dimension(NXMEMQ_,NYMEM_,NZ_), intent(out)   :: PFu
  real, dimension(NXMEM_,NYMEMQ_,NZ_), intent(out)   :: PFv
  type(ocean_grid_type),               intent(in)    :: G
  type(PressureForce_AFV_CS),          pointer       :: CS
  real, dimension(:,:),               optional, pointer     :: p_atm
  real, dimension(NXMEM_,NYMEM_,NZ_), optional, intent(out) :: pbce
  real, dimension(NXMEM_,NYMEM_),     optional, intent(out) :: eta

!    This subroutine determines the acceleration due to pressure forces, using
!  a finite volume form of the terms and analytic integrals in depth to avoid
!  the possibility of a thermobaric instability, as described by Adcroft et al.,
!  2008.
!    To work, the following fields must be set outside of the usual
!  ie to ie, je to je range before this subroutine is called:
!   h[ie+1] and h[je+1] and (if tv%eqn_of_state is set) T[ie+1], S[ie+1],
!   T[je+1], and S[je+1].
! Arguments: h - Layer thickness, in H (probably m).
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (out)     PFu - Zonal acceleration due to pressure gradients
!                  (equal to -dM/dx) in m s-2.
!  (out)     PFv - Meridional acceleration due to pressure
!                  gradients (equal to -dM/dy) in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 PressureForce_init.
!  (in)      p_atm - the pressure at the ice-ocean or atmosphere-ocean
!                    interface in Pa.
!  (out)     pbce - the baroclinic pressure anomaly in each layer
!                   due to free surface height anomalies, in m s-2.
!  (out)     eta - the free surface height used to calculate PFu and PFv, in m,
!                  with any tidal contributions or compressibility compensation.

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    rho_star    ! In-situ density divided by the derivative with depth of the
                ! corrected e times (G_Earth/Rho0).  In units of m s-2.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e ! Interface height in m.
  real, dimension(SZI_(G),SZJ_(G))  :: &
    dz, &       ! The change in geopotential thickness through a layer, m2 s-2.
    e_tidal, &  ! The bottom geopotential anomaly due to tidal forces from
                ! astronomical sources and self-attraction and loading, in m.
    dM, &       ! The barotropic adjustment to the Montgomery potential to
                ! account for a reduced gravity model, in m2 s-2.
    pa, &       ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the
                ! the interface atop a layer, in Pa.
    dpa, &      ! The change in pressure anomaly between the top and bottom
                ! of a layer, in Pa.
    Rho_cv_BL,& !   The coordinate potential density in the deepest variable
                ! density near-surface layer, in kg m-3.
    intz_dpa    ! The vertical integral in depth of the pressure anomaly less
                ! the pressure anomaly at the top of the layer, in m Pa.
  real, dimension(SZIQ_(G),SZJ_(G)) :: &
    intx_pa, &  ! The zonal integral of the pressure anomaly along the interface
                ! atop a layer, divided by the grid spacing, in Pa.
    intx_dpa    ! The change in intx_pa through a layer, in Pa.
  real, dimension(SZI_(G),SZJQ_(G)) :: &
    inty_pa, &  ! The meridional integral of the pressure anomaly along the
                ! interface atop a layer, divided by the grid spacing, in Pa.
    inty_dpa    ! The change in inty_pa through a layer, in Pa.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
    T_tmp, &    ! Temporary array of temperatures where layers that are lighter
                ! than the mixed layer have the mixed layer's properties, in C.
    S_tmp       ! Temporary array of salinities where layers that are lighter
                ! than the mixed layer have the mixed layer's properties, in psu.
  real :: rho_in_situ(SZI_(G)) ! The in situ density, in kg m-3.
  real :: p_ref(SZI_(G))     !   The pressure used to calculate the coordinate
                             ! density, in Pa (usually 2e7 Pa = 2000 dbar).
  real :: p0(SZI_(G)) ! An array of zeros to use for pressure in Pa.
  real :: h_neglect          ! A thickness that is so small it is usually lost
                             ! in roundoff and can be neglected, in m.
  real :: I_Rho0             ! 1/Rho0.
  real :: G_Rho0             ! G_Earth / Rho0 in m4 s-2 kg-1.
  real :: Rho_ref            ! The reference density in kg m-3.
  real :: Rho_anom_here      ! The in-situ density, averaged over a layer, less
                             ! Rho_ref, in kg m-3.
  real :: raL, raR           ! The in-situ density at the thickness point to the
                             ! "left" or "right" of a velocity point, averaged
                             ! over a layer, less Rho_ref, in kg m-3.
  logical :: use_p_atm       ! If true, use the atmospheric pressure.
  logical :: use_EOS    ! If true, density is calculated from T & S using an
                        ! equation of state.
  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
  real, parameter :: C1_6 = 1.0/6.0
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
  integer :: i, j, k

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  nkmb=G%nk_rho_varies
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq

  use_p_atm = .false.
  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
  use_EOS = associated(tv%eqn_of_state)
  do i=Isq,Ieq+1 ; p0(i) = 0.0 ; enddo

  if (.not.associated(CS)) call GOLD_error(FATAL, &
       "GOLD_PressureForce: Module must be initialized before it is used.")

  h_neglect = G%H_subroundoff
  I_Rho0 = 1.0/G%Rho0
  G_Rho0 = G%g_Earth/G%Rho0
  rho_ref = CS%Rho0

  if (CS%tides) then
    !   Determine the surface height anomaly for calculating self attraction
    ! and loading.  This should really be based on bottom pressure anomalies,
    ! but that is not yet implemented, and the current form is correct for
    ! barotropic tides.
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      e(i,j,1) = -1.0*G%D(i,j)
    enddo ; enddo
    do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      e(i,j,1) = e(i,j,1) + G%H_to_m*h(i,j,k)
    enddo ; enddo ; enddo
    call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp)
  endif

!    Here layer interface heights, e, are calculated.
  if (CS%tides) then
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      e(i,j,nz+1) = -1.0*G%D(i,j) - e_tidal(i,j)
    enddo ; enddo
  else
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      e(i,j,nz+1) = -1.0*G%D(i,j)
    enddo ; enddo
  endif
  do k=nz,1,-1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
    e(i,j,K) = e(i,j,K+1) + G%H_to_m*h(i,j,k)
  enddo ; enddo ; enddo


  if (use_EOS) then
! With a bulk mixed layer, replace the T & S of any layers that are
! lighter than the the buffer layer with the properties of the buffer
! layer.  These layers will be massless anyway, and it avoids any
! formal calculations with hydrostatically unstable profiles.

    if (nkmb>0) then
      tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp
      tv_tmp%eqn_of_state => tv%eqn_of_state
      do k=1,nkmb ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
      enddo ; enddo ; enddo
      do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
      do j=Jsq,Jeq+1
        call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
                        Rho_cv_BL(:,j), Isq, Ieq-Isq+2, tv%eqn_of_state)
      enddo
      do k=nkmb+1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        if (G%Rlay(k) < Rho_cv_BL(i,j)) then
          tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
        else
          tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
        endif
      enddo ; enddo ; enddo
    else
      tv_tmp%T => tv%T ; tv_tmp%S => tv%S
      tv_tmp%eqn_of_state => tv%eqn_of_state
    endif
  endif

  if (CS%GFS_scale < 1.0) then
    ! Adjust the Montgomery potential to make this a reduced gravity model.
    if (use_EOS) then
      do j=Jsq,Jeq+1
        if (use_p_atm) then
          call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), &
                                 rho_in_situ, Isq, Ieq-Isq+2, tv%eqn_of_state)
        else
          call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, &
                                 rho_in_situ, Isq, Ieq-Isq+2, tv%eqn_of_state)
        endif
        do i=Isq,Ieq+1
          dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1)
        enddo
      enddo
    else
      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * G%Rlay(1)) * e(i,j,1)
      enddo ; enddo
    endif
  endif

  ! Set the surface boundary conditions on pressure anomaly and its horizontal
  ! integrals, assuming that the surface pressure anomaly varies linearly
  ! in x and y.
  if (use_p_atm) then
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      pa(i,j) = (rho_ref*G%g_Earth)*e(i,j,1) + p_atm(i,j)
    enddo ; enddo
  else
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      pa(i,j) = (rho_ref*G%g_Earth)*e(i,j,1)
    enddo ; enddo
  endif
  do j=js,je ; do I=Isq,Ieq
    intx_pa(I,j) = 0.5*(pa(i,j) + pa(i+1,j))
  enddo ; enddo
  do J=Jsq,Jeq ; do i=is,ie
    inty_pa(i,J) = 0.5*(pa(i,j) + pa(i,j+1))
  enddo ; enddo

! Have checked that rho_0 drops out and that the 1-layer case is right. RWH.

  do k=1,nz
    ! Calculate 4 integrals through the layer that are required in the
    ! subsequent calculation.
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      dz(i,j) = G%g_Earth*G%H_to_m*h(i,j,k)
    enddo ; enddo
    if (use_EOS) then
      call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
                          rho_ref, CS%Rho0, G%g_Earth, G, tv%eqn_of_state, &
                          dpa, intz_dpa, intx_dpa, inty_dpa)
    else
      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        dpa(i,j) = (G%Rlay(k) - rho_ref)*dz(i,j)
        intz_dpa(i,j) = 0.5*(G%Rlay(k) - rho_ref)*dz(i,j)*h(i,j,k)
      enddo ; enddo
      do j=js,je ; do I=Isq,Ieq
        intx_dpa(I,j) = 0.5*(G%Rlay(k) - rho_ref) * (dz(i,j)+dz(i+1,j))
      enddo ; enddo
      do J=Jsq,Jeq ; do i=is,ie
        inty_dpa(i,J) = 0.5*(G%Rlay(k) - rho_ref) * (dz(i,j)+dz(i,j+1))
      enddo ; enddo
    endif

    do j=js,je ; do I=Isq,Ieq
      PFu(I,j,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - &
                     (pa(i+1,j)*h(i+1,j,k) + intz_dpa(i+1,j))) + &
                    ((h(i+1,j,k) - h(i,j,k)) * intx_pa(I,j) - &
                     (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa(I,j))) * &
                   ((2.0*I_Rho0*G%IDXu(I,j)) / &
                    ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
      intx_pa(I,j) = intx_pa(I,j) + intx_dpa(I,j)
    enddo ; enddo
    do J=Jsq,Jeq ; do i=is,ie
      PFv(i,J,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - &
                     (pa(i,j+1)*h(i,j+1,k) + intz_dpa(i,j+1))) + &
                    ((h(i,j+1,k) - h(i,j,k)) * inty_pa(i,J) - &
                     (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa(i,J))) * &
                   ((2.0*I_Rho0*G%IDYv(i,J)) / &
                    ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
      inty_pa(i,J) = inty_pa(i,J) + inty_dpa(i,J)
    enddo ; enddo
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      pa(i,j) = pa(i,j) + dpa(i,j)
    enddo ; enddo
  enddo


  if (CS%GFS_scale < 1.0) then
    do k=1,nz
      do j=js,je ; do I=Isq,Ieq
        PFu(I,j,k) = PFu(I,j,k) - (dM(i+1,j) - dM(i,j)) * G%IDXu(I,j)
      enddo ; enddo
      do J=Jsq,Jeq ; do i=is,ie
        PFv(i,J,k) = PFv(i,J,k) - (dM(i,j+1) - dM(i,j)) * G%IDYv(i,J)
      enddo ; enddo
    enddo
  endif

  if (present(pbce)) then
    call set_pbce_Bouss(e, tv_tmp, G, G%g_Earth, CS%Rho0, CS%GFS_scale, pbce)
  endif

  if (present(eta)) then
    if (CS%tides) then
    ! eta is the sea surface height relative to a time-invariant geoid, for
    ! comparison with what is used for eta in btstep.  See how e was calculated
    ! about 200 lines above.
      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        eta(i,j) = e(i,j,1) + e_tidal(i,j)
      enddo ; enddo
    else
      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        eta(i,j) = e(i,j,1)
      enddo ; enddo
    endif
  endif

  if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag)

end subroutine PressureForce_AFV_Bouss


subroutine PressureForce_AFV_init(Time, G, param_file, diag, CS, tides_CSp)
  type(time_type), target, intent(in)    :: Time
  type(ocean_grid_type),   intent(in)    :: G
  type(param_file_type),   intent(in)    :: param_file
  type(diag_ptrs), target, intent(inout) :: diag
  type(PressureForce_AFV_CS),  pointer       :: CS
  type(tidal_forcing_CS), optional, pointer :: tides_CSp
! Arguments: Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure containing pointers to common diagnostic fields.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module.
!  (in)      tides_CSp - a pointer to the control structure of the tide module.
  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  character(len=40)  :: mod   ! This module's name.

  if (associated(CS)) then
    call GOLD_error(WARNING, "PressureForce_init called with an associated "// &
                            "control structure.")
    return
  else ; allocate(CS) ; endif

  CS%diag => diag ; CS%Time => Time
  if (present(tides_CSp)) then
    if (associated(tides_CSp)) CS%tides_CSp => tides_CSp
  endif

  mod = "GOLD_PressureForce_AFV"
  call log_version(param_file, mod, version, tagname)
  call get_param(param_file, mod, "RHO_0", CS%Rho0, &
                 "The mean ocean density used with BOUSSINESQ true to \n"//&
                 "calculate accelerations and the mass for conservation \n"//&
                 "properties, or with BOUSSINSEQ false to convert some \n"//&
                 "parameters from vertical units of m to kg m-2.", &
                 units="kg m-3", default=1035.0)
  call get_param(param_file, mod, "TIDES", CS%tides, &
                 "If true, apply tidal momentum forcing.", default=.false.)

  if (CS%tides) then
    CS%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', G%axesh1, &
        Time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter')
  endif

  CS%GFS_scale = 1.0
  if (G%g_prime(1) /= G%g_Earth) CS%GFS_scale = G%g_prime(1) / G%g_Earth

  call log_param(param_file, mod, "GFS / G_EARTH", CS%GFS_scale)

end subroutine PressureForce_AFV_init


subroutine PressureForce_AFV_end(CS)
  type(PressureForce_AFV_CS), pointer :: CS
  if (associated(CS)) deallocate(CS)
end subroutine PressureForce_AFV_end

end module GOLD_PressureForce_AFV
